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We derive analytical formulae for the firing rate of integrate-and-fire neurons endowed with real- 
istic synaptic dynamics. In particular we include the possibility of multiple synaptic inputs as well 
as the effect of an absolute refractory period into the description. 
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I. INTRODUCTION 

In vivo neurons in cortical and other neural circuits 
experience a large background of synaptic inputs, act- 
ing as a source of noise. Noisy inputs have an impor- 
tant impact on the dynamics of neurons, making neu- 
ral responses highly variable and affecting many of their 
response characteristics 0, i, i, i, S- A fundamental 
problem is thus to determine the output statistics of the 
neuronal activity given an input noise statistics. In addi- 
tion, the knowledge of the neuronal firing properties can 
be used to explore large-scale networks using a mean- 
field approach. In this framework the stationary states 
of populations of interacting neurons are self-consistently 
obtained from their firing responses 0, 0, @] ■ This allows 
the efhcient exploration of the parameters space and the 
characterization of the various regime of functioning of 
these networks. For these reasons it appears crucial to 
have an accurate estimation of the input-ouput relation- 
ship of neurons, especially, in presence of realistic synap- 
tic currents. 

In this direction we investigate the firing rates of 
integrate-and-fire (IF) neurons. The firing frequency of 
neurons with instantaneous synaptic inputs was first ob- 
tained in Ref. [§| . The effect of synaptic dynamics has 
been studied under various assumptions [l^, [O, ■ 
Here we derive exact formulae that include additional im- 
portant features. First we take into account the presence 
of a finite refractory time, which leads to current correla- 
tions between spikes affecting the firing rate of the neu- 
ron. Indeed, neurons at their reset potential will evolve 
with a synaptic current that is still correlated with the 
positive going current that made them cross the thresh- 
old. Furthermore, we also consider the case of multiple 
inputs from different synaptic receptors, as it occurs in 
the vast majority of cortical circuits. 

The obtention of the firing rate can be recast into 



the form of a mean first passage time (MFPT) calcu- 
lation. Since the synaptic dynamics here plays a cen- 
tral role we have to consider a multi-dimensional Fokker- 
Planck equation. The nonequilibrium distribution of the 
synaptic currents at the reset potential implies that the 
proper phase space boundary conditions must be found 
self-consistently from the neuronal and synaptic dynam- 
ics. To cope with these issues we extend the rece ntly 
developed tools by Doering and coworkers [3, [H, [l6| . 
In this approach the relevant parameter is the ratio be- 
tween the synaptic current and membrane potential time 
constants. 

Similar considerations also appear in various areas of 
science when studying the escape rate from a metastable 
state. The present results are of special interest in the 
context of stochastic resonance or stochastic activation 
[TtI . [isj , where considerable attention has been paid to 
the coexistence of several colored noises with nonequilib- 
rium distributions. 



II. IF NEURONS AND SYNAPTIC DYNAMICS 

The sub-threshold neuronal dynamics of integrate-and- 
fire neurons is described by the depolarization V{t) of the 
soma, which evolves according to an evolution equation 
of the form 19] 



(1) 



with Tm the membrane time constant. The function f{V) 
governs the dynamics of the membrane voltage when no 
synaptic currents are present. For / = we have a per- 
fect integrate-and-fire neuron; for f{V) — —V we have 
a leaky integrate and fire neuron. / denotes the synap- 
tic current and g its associated input conductance. The 
synaptic current evolves according to 



(2) 
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where ^ denotes a Gaussian white noise with zero mean 
and unit variance: {£,{t)) = and (^(i)C(i')) = - 



2 



t'). The synaptic current is thus exponentially correlated 
in time with a correlation time r,: 



hm {l{t)I{t + T)) 



D 
2^ 



exp(-|r|/rs) , 



(3) 



in terms of /(t) — I{t) — fi, where fi is the average cur- 
rent, and the noise intensity D. This form of the synaptic 
dynamics arises when the neuron receives a large number 
of inputs during its characteristic time, as it is typically 
the case in vivo in the cortex [l^. In the limit of in- 
stantaneous synaptic current, —^ 0, or for times much 
larger than the correlation time Tg , the dynamics simpli- 
fies to T,ndV/dt = f{V) + fj./g + y/D/gl{t). Note that 
the adjunction of the synaptic dynamics ((2]) renders the 
voltage dynamics by itself non-Mar kovian. 

The neuron fires an action potential when the voltage 
reaches the threshold potential Vr- After an absolute 
refractory period of duration it is reset at the value 
Vr < Vt- During the refractory period no further firing 
can occur. By contrast, even when the neuron is in the 
refractory state, the synaptic current continues to evolve 
under Eq. leading to current correlations of order 
exp(— t^/ts) between spikes. This source of correlations 
must be taken into account in order to obtain a self- 
consistent input-output relationship. 

In the following we will work in the reduced variables 



z — e{I — /i) /a and v — gV , 
where we introduced the parameters 

e = \/TslTrr, and = D/2Tm 



(4) 



(5) 



Using the adimensional time i^™ = t/T,n we obtain the 
system of equations 



V = — M yv) H — z 
e 

i = -- + — 



(6a) 
(6b) 



The potential u{v) is such that u'(v) — —gf(v/g) — /z, 
where the prime denotes a derivative with respect to v. 
The threshold and reset potentials become 



gVr and rj = gVR 



(7) 



in these new variables. Henceforth we will also refer to 
the variables v and z as the voltage and the synaptic 
current, respectively. 



III. STATIONARY FIRING RATE 

The process ([6|) obeys the Fokker-Planck equation [H 



dtP 



e '^dz{z + dz) — e ^azd^ + dyu' 



P 



(8) 



for the joint probability distribution P{v,z,t). We will 
treat e as a small parameter, considering the situation 



where the synaptic time scale is smaller than the mem- 
brane time constant. This is for example the case for 
AMPA receptors, which constitute the main fast exci- 
tatory inputs in the brain and have a time constant 
■''AMPA ~ 2 ms smaller jthan the time constant of pyrami- 
dal cells, Tm ^ 25 ms [2l|. Note that the variable z keeps 
a finite variance irrespective of the value of the parameter 
e. 

Following Doering et al we have a singular pertur- 
bation problem for the quantity 



Qiv,z) 



P{v,z,t)dt. 



(9) 



The latter gives the mean time a neuron spends at points 
(u, z) before crossing the threshold. Alternatively, this 
stationary problem corresponds to the situation where we 
perform a time average over a long trajectory where the 
neuron restarts its time evolution at the reset potential 
after firing and after its absolute refractory period. After 
some transients the time spent at each point in phase 
space will reach the stationary value Q{v,z) [l^]. The 
mean first passage time before firing then reads 



+ 00 



dz 



dv Q{v, z) 



while the firing rate of the neuron is given by 

1 



(10) 



(11) 



The factor accounts for the lowering of the firing rate 
due to the time spent in the refractory state. Impor- 
tantly, we must consider the additional dependance on 
the refractory time that appears through the MFPT it- 
self. 

To obtain the mean first passage time (1101) we first gen- 
eralize the calculation of Doering et al [IJI to the case 
where the synaptic current at reset potential follows an 
arbitrary distribution pl{z). This distribution will have 
to be determined self-consistently from the neuronal dy- 
namics. The function Q{v, z) obeys the equation 



azdv + dyU Q{v, z) — —S{v — ri)fi{z) , (12) 



where we introduced the operator L = dz{z dz). Fur- 
thermore, it must satisfy the half-line absorbing bound- 
ary condition [2^ 



,z) = for z < 



eu'{e) 



(13) 



This condition stems from the observation that the prob- 
ability current J„ = {—u'(v) +az/t)Q must be positive 
at the threshold potential for firing to occur. Intuitively 
this corresponds to the fact that the potential cannot 
cross the threshold from above. 
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For further reference we here introduce the eigenfunc- 
tions of the operator L: 

Lpn{z) = -npn{z) , Pn{z) = -=^Hc„(z) , (14) 

where He„(z) = (-iyV='/2(d"/dz")e-^'/2 ^re the 
Hermite polynomials [2J]. The latter form a 
complete orthogonal basis in the inner product 
J-^ He„(z)Hem(z) exp(— z^/2) — \p2j:n\ 6nm- Hence we 
may introduce the quantities 

1 /■+°° 

c„ = -7 / fi{z)Ren{z)dz (15) 

characterizing the distribution fi. Note that cq = 
p{z) dz — 1 and that ci = l^i^)^ dz is the mean 
of the distribution p. 

We now insert the ansatz 

Q{v, z) = Qo(«, z) + eQi(w, z) + (^Q2{v., z) + . . . (16) 

into Eq. p2)) and collect terms according to powers of 
e. The zeroth order corresponds to the white noise limit 
of Eqs. (I6|), i.e., to the limit where synaptic events are 
instantaneous. We have 

LQq(v,z) = ^ (17) 

whose solution is 

QQ(y,z)=rQ{v)pQ{z), (18) 

with a yet undetermined function ro(w). At the next 
order of perturbation we find 

L(3i(i;,z) = crpi(z)9„ro(w) , (19) 

yielding 

Q\{v, z) = ri{v)po{z) - api{z)dyro{v) (20) 

with the yet unknown function ri(u). The second-order 
terms give 

oo 

LQ2(v,z) = -S{v - ri)'^CnPn{z) + pi{z)adyri{v) 

n=0 

~P2iz)a^d.uro{v) - po{z)[dy{urQ) + cr^5^ro(w)] , (21) 

where we used that zpn{z) = pn+i{z) + pn-i{z). The 
operator L is invertible only in the subspace of functions 
spanned by n > 1. Hence the term proportional to po 
must vanish, leading to the following equation for ro(w): 

a^dlro + (u'ro) = -co6{v ~ 77) (22) 

with Co = 1. Similarly, the integrability condition for 
leads to the following equation for ri{v): 

a^dln + 9„(uVi) = -ciadJiv - -q) . (23) 



The solutions of these equations that are integrable read 

ro{v) = a-2e-«(")/^' T dw' e"("')/"'e(v' - ??) 
Je 

+ A(7-2e-«(")/<T' (24) 

and 

+ cif7-ie["(")-"(")l/"'e(?7-w), (25) 

where 8 is the Heaviside function. Remarkably the first- 
order correction to the MFPT, and hence to the firing 
rate, only depends on ci, that is on the mean of the 
current distribution p{z). Note that the MFPT can 

be expressed as (T) = J^^[ro{v) + eri{v) + ■ ■ ■]dv as 
the functions rn{v) are the reduced densities, r„(w) = 

J^^ Qn{v,z)dz. 

We now have to determine the factor ci and the con- 
stants A and B self-consistently. In this regard, we first 
consider the phase space boundary condition p3p . The 
latter was solved in Ref. and can be summarized by 
the condition 

r{e) = eaar'{6) (26) 

where a = — C(l/2) ~ 1.46 . . . with ( the Riemann zeta 
function. The crucial point to observe is that this bound- 
ary condition remains valid in presence of an arbitrary 
current distribution p since, as shown by Eq. (j25|) . the 
presence of a non-vanishing mean current, ci 7^ 0, only 
affects the residence time in the region v < rj, leaving 
the region rj < v < 9 unchanged. Hence the constants 
A and B take the values yl = and B ^ a. The factor 
a is referred to as the Milne extrapolation length as it 
characterizes the non-vanishing value of the probability 
distribution at the threshold potential. 

We now have to evaluate the value of the mean current 
ci at the reset potential. It must be determined self- 
consistently as the synaptic current remains correlated 
between spikes due to its finite time constant. Accord- 
ingly, the probability distributions at the threshold and 
reset potentials are related as follows: 

/+00 
Giz,z',Tr)Qi9,z')dz' (27) 
-00 

where G{z, z', t) is the Green function of the synaptic 
current, which is given by a normal distribution for z of 
mean z'exp(— i/e^) and variance 1 — exp(— 2t/e^). Con- 
sequently we have that 

/+00 
z'Q{e,z')dz' , (28) 
-00 

where we used Eq. ^ with 00(2, z', Tr)dz = 

z' cxp(— t^/ts). The mean reset current is thus related to 
the mean threshold current. Again, we can verify that. 
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at first-order in e, the expression for the probabihty flux 
given in Ref. [15] remains unchanged in the present sit- 
uation, as aheady suggested by Eq. ([25|) . Hence we find 



ci = aexp(-Tr/rs) 



(29) 



at first-order in e. When the refractory time is large com- 
pared to the synaptic time constant the current will relax 
to its stationary distribution ii{z) = /^/\/2tt so that 
c„ = for n > 1. When the refractory time vanishes cor- 
relations are maximal since the current distribution fi at 
the reset potential exactly matches the current distribu- 
tion at the threshold, resulting in a contribution ci — a. 
The same method can be applied to obtain the full dis- 
tribution of the synaptic current. 

Combining these results we obtain the firing rate pT|) 
of IF neurons. An important case is the leaky integrate- 
and fire neuron for which the mean first passage time 
takes the compact form 



(T) 



e'" [erf(w) + l]dw 



(30) 



where rj* = (Vr - fj,/g)/V2a + eaexp{-Tr/Ts)/V2, 0* = 
(Vt — [ilg)l\f2a + ea/-\/2, and erf is the error function 
[24| . This formula interpolates between the two-limiting 
cases of zero and infinite refraction time previously con- 
sidered in Refs ll] and pj|, respectively. Similar expres- 
sions can be obtained for nonlinear IF neurons as well. 

The firing rate pT|) as a function of the refractory pe- 
riod is depicted in Fig. [1] along with the result of Monte- 
Carlo simulations. The firing rate is compared to the 
situation where the neuron is unresponsive during the 
refractory period but the corresponding current correla- 
tions are neglected [i.e., formula (fTTj) but with (T)^^]. 
The corrections due to the effect of the refractory state 
on the MFPT are significant, even at moderate firing 
rates. 

To investigate the validity of the expansion (|16p . we 
numerically obtained the firing rate as a function of the 
ratio Ts/Tm, which is plotted in Fig. [21 The agreement 
with the first-order formula (j30p holds for the region 
Ts/Tm < 0.1. For ratios Ts/Tm up to 1, we found that 
the firing rate could be well fitted by considering the ef- 
fective threshold 



(I/i^-M/.g + 0.1375e2) 



0.225e^ (31) 



V2(T V2 

as seen in Fig. [51 The effective reset potential i]* appears 
unaffected by the second order corrections, as observed in 
Ref. [10] . In particular, the dependance on the refractory 
period remains identical. 

IV. MULTIPLE SYNAPTIC INPUTS 

We further consider the case where several synaptic 
inputs are present, i.e.. 




T (ms) 

FIG. 1: Firing rate of leaky integrate-and-fire neurons as a 
function of the absolute refractory period. The dashed lines 
represent the approximation of stationary synaptic current 
(see text). Dots depict the result of Monte-Carlo simulations. 
The reset and threshold potentials take the values Vt = 15 
mV and Vr — 20 mV, respectively (the resting potential is 
mV). The membrane time constant Tm = 25 ms and the 
synaptic time constant Ts = 2 ms so that e ^ 0.283. The 
mean potential n/g — 15.5 mV and its standard deviation 
a — 5 mV for the lower curve (black), and n/g — 16.5 mV 
and a — 6 mV for the upper curve (red). 
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FIG. 2: Firing rate of leaky integrate-and-fire neurons as a 
function of the ratio Ts/rm- Monte-Carlo simulations (sym- 
bols) are compared to formula (lll|l with the effective thresh- 
old (|3ip and Tr = (dashed lines). Different curves corre- 
spond to different means and variances of the synaptic input. 



where the currents Ik obey evolution equations of the 
form ([2]) albeit with different time constants r^^ , means 
We here assume that the synap- 



with c of order 



I = h+l2 



(32) 



/ife, and variances a'j. 
tic time constants satisfy Ts^/ts, 
unity, c = 0(1), so that Ts^. are of order e: -^/tVi/tv^ = 
c^Ts^jTm = e. The calculations of the previous section 
can be extended to this case by considering the three- 
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dimensional phase space Fokker-Planck equation. Insert- 
ing the eigenfunctions Pnizi)Pmiz2) into the expansion in 
powers of e, we obtain that the boundary condition now 
reads €C7\r'{0) = r(0), where cr^ — al + a^- We could 
not obtain the Milne extrapolation length analytically, 
which was determined numerically to be A ~ 1.4. It re- 
places the factor a = — C(l/2) ~ 1.46 ... in the present 
situation. 

The white noise limit is obtained from Eq. ((24| with 
C = and = af + and where the potential u 
incorporates both input currents, p — fJ-i + ^2- The first- 
order correction reads 

=Ae["(°)-"("^l/"', (33) 

with the corresponding correction to the MFPT given 
by e ri{v)dv (for simplicity we did not include the 
refractory state in this analysis). Remarkably, this first- 
order correction only depends on the total noise strength 
and not on the ratio of the synaptic time constants. 
Accordingly, as regards the firing rate, we may treat dif- 
ferent synaptic currents as an effective single input with 
parameters p — p.i + p,2 and — + a"^, even when 
they present different time courses. This effect is valid 
at first-order in the parameter e; the case of very long 
synaptic time constants compared to the membrane inte- 
gration time was studied in Ref. [12]. The generalization 
to ri > 2 synaptic inputs is straightforward and does not 
change this result. 



The finiteness of the refractory period of the neuron 
provides a new source of correlations interacting with the 
synaptic dynamics. Since the synaptic dynamics has a fi- 
nite relaxation time, neurons at the reset potential may 
still be correlated with the positive current that made 
them cross the firing threshold, resulting in an increased 
firing rate. This mechanism is important for fast excita- 
tory and inhibitory synapses such as AMPA and GABA^i 
synapses as their time constants is in the range of 2 to 5 
ms [2l|, comparable to the absolute refractory period of 
cortical neurons, 1 — 5 ms. 

The presence of multiple synaptic inputs has been con- 
sidered as well. In this case the Milne extrapolation 
length must be determined numerically to complete the 
theoretical analysis. Importantly the firing rate only de- 
pends on the variances of the synaptic inputs and not on 
their relative time constants. This allows us to consider 
the combined effect of excitatory and inhibitory synaptic 
inputs from, say, AMPA and GKQAa receptors. 

These analytical results are crucial to assess the effect 
of noise on neuronal dynamics. For instance they provide 
accurate expressions that can be used in the mean-field 
exploration of large-scale neuronal networks @, 0, HI ■ In 
particular, the network collective properties, such as the 
stability of the low activity spontaneous state or of the 
persistent "memory" states, will depend on the transfer 
function. More generally, the study of nonequilibrium, 
multiple colored noises appears in a wide range of situa- 
tions of concern, from chemical networks to biology. 



V. CONCLUSIONS 

We have obtained analytical expressions for the sta- 
tionary firing rate of integrate-and-fire neurons endowed 
with realistic synaptic dynamics. Precisely we have in- 
cluded two salient features into the description: the pres- 
ence of a finite refraction time on the one hand, and the 
presence of multiple synaptic inputs on the other hand. 
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